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Statistical Challenges in Microrheology *^ 
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Abstract 



Microrheology is the study of the properties of a complex fluid through the diffusion dy- 
namics of small particles, typically latex beads, moving through that material. Currently, it is 
, the dominant technique in the study of the physical properties of biological fluids, of the ma- 

terial properties of membranes or the cytoplasm of cells, or of the entire cell. The theoretical 
underpinning of microrheology was given in Mason & Weitz (Physical Review Letters; 1995), 
who introduced a framework for the use of path data of diffusing particles to infer viscoelastic 
properties of its fluid environment. The mult i- particle tracking techniques that were subse- 
r-fl ■ quently developed have presented numerous challenges for experimentalists and theoreticians. 

This paper describes some specific challenges that await the attention of statisticians and 
applied probabilists. We describe relevant aspects of the physical theory, current inferential 
efforts and simulation aspects of a central model for the dynamics of nano-scalc particles in 
viscoelastic fluids, the generalized Langcvin equation. 



(N ■ 1 Introduction 

> 

\ Broadly speaking, rheology is the study of the flow and deformation of soft matter. Classically, 
OS " a primary goal of rheologists is to develop constitutive laws for colloids, gels and non-Newtonian 
fluids. The field is filled with amusing household examples of mathematically exotic behavior: 
shear-thinning (for example, the way one can help ketchup pour out of a bottle by shearing it 
with a knife), rod climbing (the way cake batter rises up a stirrer as it is mixed), and shear- 
thickening (the way some fluids, such as the children's toy Oobleck, exhibit properties of a solid 
when subjected to shear). Many biological fluids, such as blood, mucus and the cytoplasm of 
cells, exhibit non-Newtonian properties. Models of diffusion in biological fluids have been devel- 
oped for pharmaceutical [651 E3 H2] and medical [21] use. These applications exploit both local 
properties, such as polymer mesh size [H] and viscoelasticity [301 [26], an d global properties, such 
as heterogeneity of both viscosity and layer thickness |41] , 



*The first author was supported in part by the Louisiana Board of Regents award LEQSF(2008-ll)-RD-A-23. 
The third author was supported by the NSF grant DMS-1100281, the Cystic Fibrosis Foundation (HILL08I0), and 
the NIH (R01-HL077546-01A2). The fourth author was supported in part by the NSF grant DMS-0714939. All of 
the authors would like to thank Gregory Forest. 

^ AMS Subject classification. Primary: 60G15, 82B31, 62P10. Secondary: 62M09, 42C40. 

^Keywords and phrases: Gaussian processes, microrheology, non- Newtonian fluids, generalized Langevin equa- 
tion, local Whittle, fractional processes, simulation, wavelets. 
§ Mathematics Department, Tulane University. 
1 Mathematics Department, University of Florida. 

" Cystic Fibrosis Pulmonary Research and Treatment Center, University of North Carolina at Chapel Hill. 
"Department of Statistics, Pennsylvania State University. 

^Correspondence to: John Fricks, Department of Statistics, 326 Thomas Building, Penn State University, Uni- 
versity Park, PA 16802-2111. E-mail: fricks@stat.psu.edu. 



1 



The rheological study of such biological materials poses a unique practical problem though. 
Whereas classical rheometers can probe materials by applying both local and global shears and 
stresses, biological materials tend not to be available in sufficient quantities to apply traditional 
tests. (That is to say, experimenters find it very difficult to come across liters of mucus on 
demand.) To confront this constraint, techniques have been developed to probe media locally 
by embedding and tracking manufactured microparticles [M]. Dragging a magnetized microbead 
through a medium is considered active microrheology, while simply placing a bead and recording 
its diffusive motion is considered passive. For biological materials, the latter is often preferred 
because the active process may destroy the medium as it is being measured. 

At first blush, it seems unlikely that passive observation of diffusing particles can reveal much 
about a fluid's microstructure. However, an important step in this direction appeared in the 
1995 work Optical Measurements of Frequency-Dependent Linear Viscoelastic Moduli of Complex 
Fluids, by Mason & Weitz (MW). In that work, the authors introduced a theoretical framework 
for how to use microparticle path data to infer viscoelastic properties of its fluid environment. 
The multi-particle tracking techniques that were subsequently developed have presented numerous 
challenges for experimentalists and theoreticians. In particular, the path of a single particle fol- 
lows a stochastic process; this contrasts traditional rheological techniques to study non-Newtonian 
fluids, which are essentially deterministic phenomena. Our goal is to indicate some specific chal- 
lenges arising from these stochastically driven experiments that await the attention of statisticians 
and applied probabilists. We will focus in particular on the background and current practices for 
inferring what is known as a viscoelastic fluid's complex viscosity. Where possible we will also 
indicate long-term challenges that exist for other local properties, such as directly characterizing 
a fluid's microstructure, as well as global properties like quantifying spatial heterogeneity. 

In Section 2, we summarize the MW approach. The method relies heavily on the theoretical 
mean-squared displacement (MSD), E[A 2 (i)], of embedded particles. Figure [fl from [65], is 
representative of an important group of papers [IU [48j |45j [67] that rely on MSD concepts. The 
authors display an ensemble of what are called pathwise or particle MSDs, which are notably 
nonlinear. The pathwise MSD is essentially a method of moments estimator for the theoretical 
MSD (see equation (|18l) for a definition). The phenomena of non- linear MSD is referred to as 
anomalous diffusion. There are variety of theoretical hypotheses for what can cause anomalous 
diffusion. The most prominent explanations include the presence of obstacles, caging, binding, 
memory \62\ [6"T| WQ . All of these postulates result in a particle's behavior having memory, and 
so the MW approach, which coarse-grains these microstructure concerns, is quite general. They 
use the generalized Langevin equation (GLE) as a universal model. 

In the present work, we will follow MW and use the GLE as a starting point. In Section [21 
we show conditions under which the GLE supports both diffusive and subdiffusive behavior. In 
addition, we will discuss the physical motivations for the experimental data that is being collected 
and how this data should connect to specific stochastic models. In Section [3] we note that the 
pathwise MSD techniques are subject to substantial variability, which is not well-characterized. 
This approach may not be the most efficient method for extracting information about the un- 
derlying model. This is presently overcome by ad hoc methods, but in Section [3] we discuss the 
advantages of using spectral methods such as a local Whittle estimator. First, the estimator pro- 
vides a useful rigorous test for negating a null hypothesis that a particle path is simply Brownian 
motion or other purely diffusive process. Second, since the local Whittle estimator is both more 
efficient and robust for detecting subdiffusivity, it can be used to quantify heterogeneity among 
particle paths. Finally, in Section HJ we address simulation methods including a presentation of a 
wavelet-based algorithm for the generation of fractional Brownian motion driven GLE. 
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Figure 1: At left, a snapshot from Dawson et al [16] featuring "pathwise MSDs" for particles 
of sizes 100 am, 200 nm and 500 nm respectively diffusing in sputum from a human with cystic 
fibrosis. In the center and right are figures from Valentine et al [67] displaying the ensemble 
MSD (center) and a collection of pathwise MSDs (right) of microparticles in glycerol, agarose and 
F-actin, respectively. These data are presented by microrheologists to demonstrate that complex 
fluids are viscoelastic (mucus, agarose and F-actin) and heterogeneous (mucus and agarose). 
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2 Physical Models 



2.1 Passive microrheology and Newtonian fluids 

For Newtonian fluids such as water, the rheological property of interest is the viscosity n, and it 
is a classical observation that the viscosity of a fluid can be inferred by measuring the MSD of 
suspended diffusing particles. In this section, we will develop this argument with an eye towards 
the follow-up Section 12.21 m which we describe how these notions generalize to non-Newtonian 
fluids. 

Whereas mathematicians and statisticians tend to think of Brownian motion in terms of its 
statistical properties - stationary, independent and normally distributed increments with variance 
proportional to the length of the time increment - physicists and engineers often employ a New- 
tonian framework. The construction of the celebrated Langevin equation relies on the balance of 
forces (F = ma): 

miX — -Ffriction -^external -^thermal) (1) 

where X is the position of the particle and m is its mass. For a spherical particle in a viscous 
fluid, the force due to friction is given by the Stokes drag law, Ff r i ct i on = GnrnX, where r is the 
radius of the particle and n is the viscosity of the fluid. The external force is often expressed as 
the gradient of some potential function, internal = — V$(X). In a laboratory setting, one might 
use a quadratic potential to account for the influence of an optical trap [TT], for example. In 
polymer physics and molecular dynamics simulations, which involve many particles, it is assumed 
that there is an intrinsic configuration potential that encodes the energy due to bond lengths and 
bond angles among the particles in the system (60] . 

The remaining thermal force term is meant to summarize the influence of the successive 
impacts of fluid molecules on the submerged particle. The impacts are assumed to be independent. 
Because they are frequent it is expected that a central limit theorem will hold and therefore 
the sum of their effects should be Gaussian. Furthermore, the variance of the sum should be 
proportional to both the temperature T of the system and the viscosity of the fluid rj, which both 
govern the velocity of the fluid molecules and the frequency of impact with the foreign particle. 
This observation that the variance of the fluctuations and the magnitude of the dissipative force 
of friction should should be proportional to each other, because they are both caused by the 
same physical entities in the fluid molecules, dates back to Einstein and is called the Fluctuation- 
Dissipation Theorem. 

In modern mathematical notation, physical Brownian motion X(t) in one dimension is written 
as the integral of the velocity process V(t), which is itself the solution to a stochastic differential 
equation (SDE), 

ix(t) = V(t), mdV(t) = ->yV(t)dt - V$(X(t))dt + yj2k B T^dW(t), (2) 

where k B is Boltzmann's constant. Often the SDE for the velocity process is known itself as the 
Langevin equation. 

For this section, we will restrict our attention to the free diffusion case, when §{x) = 0, 
yielding an Ornstein-Uhlenbeck process for V{t). Also, we will make calculations for diffusion in 
one dimension, but the calculations easily generalize. We take for the initial condition X(0) = 
and V(0) = v where the random variable v is drawn from the stationary distribution of V(t), 
N(0, k B T '/7m). One can write an exact solution of the free velocity equation in Duhamel form, 

V(t) = e-^v + f 'e-^dW(s), 
™ Jo 
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but more important is the autocorrelation function p(t) := E \V(to + t)V(io)]. Since the initial 
condition is chosen from the stationary distribution, the auto-correlation does not depend on to 
and is given by 

p{t) 



ksT _j_ t 



m 

Meanwhile the position process satisfies E [A(i)] = 0. Furthermore, employing Fubini's the- 
orem and linearity of expectation, we have a relationship between the MSD and the velocity 
autocorrelation function 



r rt r-t 



V{t')dt' / V(s')ds' 



Jo 



I' 



E [X 2 (t)] = E 

This simplifies to 

E [X\t)} = 2 J^Lf™(l- e -W)dt> = 2 J^(t+™(l- e -&)) 
m J 7 7 V 7 J 



p(t' - s')dt'ds'. (3) 



(4) 



If {X n (t)}t>o, n G {1, . . . , N}, are a collection of independent physical Brownian motions as 
would be the case for a Newtonian fluid assuming the particles were sufficiently spaced as to not 
interact with one another, then it follows that the empirical ensemble-averaged MSD 



(X 2 (t)):=j-J2(X n (t)-X n (0)) 2 (5) 

n=l 

will satisfy 

= m:, (6) 

t->oo t Sirrrj 

where we have written the drag 7 in terms of the viscosity of the fluid r]. Because one can control 
for the particle radius r by independent means in a laboratory setting, particle tracking can give 
a robust means to measure a fluid's viscosity. The equation © is known as the Stokes-Einstein 
Relationship and is crucial to the physical basis for passive microrheology. 



2.2 Passive microrheology and non-Newtonian fluids 

A hallmark property of non-Newtonian fluids is a nonlinear response to applied stress. A shear- 
thickening material like Oobleck, corn starch suspended in water, will exhibit a super-linear 
response to increasing applied shear force, while shear-thinning materials, like ketchup or peanut 
butter, will initially resist like a solid, but will eventually dramatically yield when sufficient force 
is applied. 

A standard procedure for capturing this nonlinear response is to apply a very small oscillatory 
shear stress. The applied force must be small because dramatic deformation events cannot yet be 
rigorously characterized. However a sufficiently small applied sinusoidal stress with magnitude 70 
and angular frequency co will induce an oscillatory response of the form [33] 

0(t) = 7o[G'(w) sin(W) + G"(cj) cos(^)]. (7) 

The functions G' and G" are respectively known as the storage and loss moduli. These are so- 
named because they capture the elastic and viscous components of the material in the following 
sense. If a material is purely elastic, the induced stress will oscillate perfectly in phase with the 
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Figure 2: A cartoon depiction of the response by viscous, elastic and viscoelastic fluids to os- 
cillatory stress. The elastic response is in phase with the stimulus, while the viscous response 
is out of phase by ir/2 radians. The viscoelastic response is characterized by a superposition of 
phase-shifted sinusoidal curves. The respective magnitudes of the elastic and viscous responses 
are often expressed as the storage and loss moduli G' and G" as displayed in Equation ([7]) and 
Figure [3j 



applied force. On the other hand, a purely viscous fluid will oscillate tt/2 radians out of phase. 
See Figure [2] for a visual representation of this phenomenon. 

The complex function G*(cj) = G'(lj) + iG"(uj) is referred to as the complex shear modulus 
and it is directly related to the natural generalization of the Newtonian viscosity, namely, the 
complex viscosity rf(uj). This relationship is given b}Q [64] 

G*(w) = iurfiuj). 

It is at this point that we rejoin the story of microrheology. 

The fundamental insight provided by MW is that the natural Generalized Stokes-Einstein 
Relation (GSER) should be written in the Laplace domain. Denoting the Laplace transform 
by Jzf[/(-)](z) = f(z) := J °° e~ zt f(t)dt and solving for the behavior of the MSD, equation ([6]) 
becomes Jl?{(X 2 )}(z) = ksT '/Ztttz 2 ^. Therefore the hypothesized form for the GSER is [47] 

*{<**>> < 8 > 

In the next section, we will offer a slightly revised version of MW's development for this form 
of the GSER. It follows from postulating that the Generalized Langevin Equation (see next 
section), first proposed by Mori and Zwanzig for the motion of particles in a heat bath |70j . is 
an appropriate model for diffusion in a viscoelastic environment. At the heart of this argument 
is a connection between the complex viscosity and an assumed memory kernel experienced by 
diffusing microparticles. 

Before moving on, though, we wish to emphasize a detail about the use of the approximation 
notation in (|8|) rather than equality (compare to similar statements in [47] and |64] ) . Most 
importantly, mass has been neglected. From this, two issues arise that have a direct impact 
on the development of statistical methods. First, the Stokes-Einstein relationship is itself an 
asymptotic relationship. One can view it as the result of taking the mass to zero or as waiting 
a sufficiently long time for inertial effects to become negligible. These asymptotic limits are 

1 lt is an idiom of the rheology community to use an asterisk to denote a Fourier transformed function. Because 
the notation is so entrenched we use G* and 77* but for other instances of the Fourier transform we use a circumflex, 
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Figure 3: Complex storage and loss moduli reported in the original paper by MW |47| . When 
G' is larger than G" , the material acts more like a solid than a fluid. When the sizes are reversed 
the material is more fluid-like. For complex fluids, whether the material is more solid or fluid 
can change depending on the frequency of applied stress. Note that the behavior of these moduli 
can either approach zero or approach a nonzero constant as the frequency u) tends to zero. This 
rheological characterization can be inferred from the asymptotic MSD of embedded microparticles. 
This connection is laid out in Section 3, equations (fT7|) . (f22l) and (f23j) . 



equivalent for viscous diffusion, but they are not necessarily equivalent for viscoelastic diffusion. 
As was demonstrated in |49j . the zero-mass limit may be singular, creating a term that can be 
lead to a misinterpretation of any observed transient anomalous diffusion. Second, as we will show 
in the next section, any large-time limit will neglect important transient features in the data. An 
entire family of classical non-Newtonian fluids known as Maxwell fluids will appear identical in a 
large-time MSD limit. 

2.3 The Generalized Langevin Equation 

As is pointed out by Kubo [39], embedded in the development of the Langevin equation ([2]) is 
an implicit assumption that there is time scale separation between the dynamics of the large 
foreign particle and those of molecules comprising the surrounding medium. However, in the 
viscoelastic regime, one cannot assume that the influences of the various fluid molecule impacts 
are independent of recent activity by the foreign particle. This manifests itself in two ways: first, 
this introduces a non-trivial autocorrelation structure in -F t h e rmai; second, the friction coefficient 
should no longer be constant. This non-Markovian property is modeled by way of a memory 
kernel, which we denote by T(t), and based on physical considerations, the fluctuation-dissipation 
theorem implies that the memory kernel should be the same for both thermal fluctuation and 
friction. The generalized Langevin equation (GLE) for the velocity process is therefore written 

m —V(t) = -i I T(t- s)V(s)ds-V$(X(t)) + F(t), (9) 
at J —oo 

where F(t) is a stationary Gaussian process with autocorrelation function 

E[F(t)F(s)] = k B T 1 T(\t-s\). (10) 

The lower limit of integration is — oo to emphasize that the velocity process is stationary. In this 
section we will continue to take the potential &(x) = for expository purposes, and we will only 
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consider 1-d dynamics. It is generally assumed that the velocity processes corresponding to each 
coordinate of a diffusing particle are independent. 

We redefine the drag parameter to be 7 = 6irrr]Q, where r/o is the viscosity of the fluid when 
there is no applied stress. The crucial physical assumption then is that the complex viscosity i]* 
is related to the Fourier transform of the memory kernel via the relation T(uj) = rj*{uS)/r)Q. For 
example, because the complex viscosity is constant for purely viscous fluids, the memory kernel 
should be a Dirac ^-distribution. Indeed, the GLE reduces to the Langevin Equation for this 
choice of V. 

For the simplest viscoelastic fluids, known as Maxwell fluids, V has the form of a Prony series: 

N 

r(t) = J2 c - e ~ Xnt - ( u ) 

n=l 

As an example, in the original paper of MW [37], the authors introduce a proof of principle 
for suspended silica particles in ethylene glycol. They effectively assumed prior knowledge of 
four distinct relaxation times (which have the form A" 1 ) and used the Laplace transform of the 
empirical MSD to fit the coefficients of the terms. 

When the number of modes is extremely large, in a polymer melt for example, such a Prony 
series can mimic a power law [361 09]. I n this case, it is customary to assume that the thermal fluc- 
tuations are given by fractional Gaussian noise (fGn). The latter is the increment process of frac- 
tional Brownian motion (fBm), which is the unique Gaussian, self-similar, stationary-increment 
process (see [23]). The associated memory kernel has the form 

T(t) = 2H(2H -l)\t\ 2H ~ 2 , t^O, -<H<1, (12) 

where H, called the Hurst parameter, characterizes the fGn/fBm. The kernel (|12p gives rise to a 
qualitatively different regime in terms of the behavior of suspended microparticles, with important 
consequences for the inferred complex viscosity. See Sections 13.31 and 14.21 

In general, because T appears as the covariance of a stochastic process, it is required to be 
a positive definite function. By Bochner's theorem, this means that there exists a real-valued 
positive Borel measure T such that T(t) = J R e~ ltu> T(du)). Deterministic linear integrodiffer- 
ential equations of a similar type to ([9]) have been well-studied [28]. If the thermal fluctuation 
is sufficiently reg for example in C(M.) or Lj oc , then Q has a path-by-path interpretation. 

Such solutions can be written in variation of constants form 

V(t) = f X (t-s)F(s)ds, (13) 

where x(i) is known as the differential resolvent in the integral equations literature [28] and as 
the susceptibility in the physics literature [2]. This function x is the fundamental solution to the 
homogeneous part of ([9]) in the sense that it satisfies 

d /■* 

m-x(i) = -7 J T(i - s)x(s) + S(t) 

2 One can still make sense of the GLE when the thermal fluctuation F is rougher. For example, when the forcing 
term is a fractional Gaussian noise there is a singularity at in the memory kernel that results in F only being 
well-defined in the sense of distributional derivatives. In Section [4] we briefly recount the approach proposed by 
Kou |36| where the spirit of the following discussion still holds. 
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Taking the Fourier transform of this equation yields the following frequency-space representation 
for x, 



To complete our characterization of the velocity process, we note that since the thermal fluctua- 
tion term is Gaussian and stationary, the linear transformation (|13p will yield another Gaussian 
stationary process. A calculation then reveals that p(t — s) = IE [V(i)V(s)] satisfies 

\miuj + 7I (u!)\ 

This sets the fundamental relationship between the velocity autocorrelation function and the 
complex viscosity. 

In order to derive the GSER ([8]) we observe that the mass of these microparticles is vanishingly 
small. Setting the mass equal to zero in the Laplace domain is equivalent to taking a certain 
weak limit of a sequence of solutions to the GLE with decreasing parameter m that may neglect 
singular correction terms [39] . If we nevertheless permit ourselves to take this step for expository 
purposes, (fLij) reduces to p(co) = ksT /^T(uj) = kBTrjo/jri*^). On the other hand, the sequence 
of equalities (|3j) implies that the Laplace transform of the autocovariance is p(z) = z 2 J2?{(X 2 )}. 
Combining these observations yields the GSER (jSJ), which expresses complex viscosity in terms 
of the MSD. 

2.4 The Generalized Langevin Equation and Anomalous Diffusion 

One of the most frequently reported properties in single particle tracking experiments is that the 
MSD scales sublinearly over time. In fact, this is the major qualitative difference between what 
one expects to see from diffusion in a Maxwell fluid versus diffusion in a polymer melt. 
To be precise, we say that a stochastic process X(t) is 

subdiffusive I E\X 2 (t)} f % 

asymptotically \ diffusive > if lim — = < a } , (15) 

superdiffusive J [00 

where a, when finite, is the diffusion coefficient. 

Morgado et al [51] demonstrated that solutions to the GLE can exhibit any of the above 
behavior depending on the Laplace transform of the memory kernel. Suppose for example that 
X(t) is a solution to ([9]) with a memory kernel T that satisfies T(s)/s 2 ds < 00 for any z > 0. 
One can show that if 

{00 I ( subdiffusive \ 

a 2 > , then the process X(t) is asymptotically < diffusive > . (16) 
J [ superdiffusive J 

To see this, recall the relationship between the MSD of the position process and the ACF of the 
velocity process ([3]): 

t' 



E[X 2 {t)}=2[ f p{t',s')dt'ds', 
Jo Jo 
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where p(t, s) = E [V(i)V(s)]. By the final value theorem, it follows from ([3]) that 



lim 



lim z££ 



E [^ 2 W] 




t— >oo 



t 



z->0 



t 



The conclusion (116j) then follows from (114p . 

This proposition implies that diffusion in Maxwell fluids is asymptotically diffusive, while 
diffusion in power law fluids will be asymptotically subdiffusive. This dichotomy plays a large 
role in the discussion of inference in the next section. In the subdiffusive case, see Section 13,31 an 
important challenge is in identifying the asymptotic subdiffusive exponent, which in turn dictates 
the near zero behavior of the complex viscosity. In the diffusive case, see Section 13.21 one must 
capture transient phenomena in order to distinguish among various model fluids. 

2.5 Modeling Challenges 

The preceding discussion has been limited to modeling ensembles of independent particles diffus- 
ing in a highly idealized version of biological fluids. We have neglected tremendously important 
factors, including surface chemistry of the particles, rigid obstacles in the medium, interaction 
with boundaries, and heterogeneity in general. For an example of the impact of such factors, we 
mention recent work by the Hanes lab, who exploit surface chemistry effects for pharmaceutical 
applications. By coating particles in polyethylene glycol (PEG), the particles are rendered mu- 
coinert and are observed to move much more quickly |65[ [42] . Interestingly, the method reveals 
certain details concerning the size of gaps in the polymer-mesh that constitutes human mucus. 
Because coated 200 nm beads diffuse faster than similarly coated 500 nm beads, this provides 
evidence that the effective mesh size is somewhere on this scale. 

The Hanes lab discovery demonstrates an important principle of passive microrheology. All 
observations are properties not of the fluid medium, but of this particle in that fluid. Therefore, 
when trying to capture rheological properties of the fluid itself, one must aggregate observations 
from multiple experimental techniques. This is one reason why two-point (and n-point in gen- 
eral) rheology is vital to the future of the field. Multiparticle tracking has the capacity to yield 
information about both local and global properties of a viscoelastic fluid. At a local level, it is 
hypothesized that fluid-particle dynamics can vary greatly with distance. For example, in cystic 
fibrosis mucus studies, it is hypothesized that depending on the surface chemistry of the particle, 
there may be a depletion layer of the mucosal network in the immediate vicinity of the particle. 
This may have a profound effect on whether a substantial number of individuals in an ensem- 
ble become immobilized due to binding. Determining this local environment is impossible with 
single bead observations though, so several recent efforts have sought to characterize the cross- 
correlations between beads |44 [ 132 1 [31] . Furthermore, depending on how rigid the fluid structure 
is, particles may "talk" to each other across long distances, as is the case for particles suspended 
in gels. 

Understanding long-range interactions among particles reveals a final inadequacy in the cur- 
rent theory. Still missing is a comprehensive model for the interaction between a viscoelastic 
medium and submerged microparticles that has the same level of sophistication that the stochas- 
tic immersed boundary method [4] provides for viscous diffusion. A robust model, coupled with 
accurate and efficient simulation techniques, are a prerequisite for characterizing whether or not 
variation seen in a given particle ensemble is greater than what should be expected from a ho- 
mogeneous medium. We give an example in the next section for how difficult this has proved for 
simple viscous fluids, revealing how much is left to be done. 



10 



3 Statistical Inference 



3.1 Pathwise Mean Squared Displacement 

The data displayed in Figure Q] is representative of numerous articles in rheological journals. In 
displaying an ensemble of pathwise MSDs (see equation (|18p ). the authors are communicating 
two important features simultaneously: the behavior of a thermally fluctuating microparticle in a 
viscoelastic medium is plainly distinct from classical diffusion, and furthermore the fluid medium 
itself can be highly heterogeneous. 

Qualitative evidence for heterogeneity can be found in the high variability of the MSD plots 
of repeated experiments. If one could propose a credible model for viscoelastic diffusion, then 
the empirical spread of pathwise MSDs can be compared to a simulated distribution. When the 
former is significantly larger than the later, this is evidence for fluid heterogeneity. An early 
attempt in this direction was performed by Valentine et al |67j . In this work, the authors used 
higher order statistics to compare empirical pathwise MSD curve distributions of microparticles 
in glycerol, agarose and F actin to the spread generated by simulated Brownian motions. 

The desire to see heterogeneity compels the use of pathwise statistics to distinguish microrhe- 
ological paths from classical diffusion. The approach is successful in that often the pathwise MSD 
curves are clearly sublinear and sometimes appear to go through a series of transitions. One 
often observes an initial inertial regime where there is ballistic (slope on a log-log plot is two) 
growth |25l 149). Over an intermediate time scale subdiffusive (slope less than one) behavior is 
observed, followed by large time scale behavior that may be subdiffusive, diffusive, or superdiffu- 
sive. 

As discussed in the previous section, the subdiffusive and diffusive long-time asymptotic be- 
haviors are natural outcomes for the GLE with power law and exponential tail memory kernels, 
respectively. While it is true that the GLE also supports super-diffusive behavior, our preliminary 
analysis of experimental data indicates that ballistic large time behavior results when there is a 
dominant drift that has not been correctly removed in the statistical analysis. To this end, it is 
also possible to introduce artificial subdiffusive behavior by removing a mean inappropriately, and 
so great care must be taken when interpreting the large time regime. 

For the purposes of the present discussion, we will assume that in a given time window /, we 
can express what one might call the "local" MSD in the following form, 



The microparticle is said to be diffusive, subdiffusive or superdiffusive if the MSD exponent a is 
equal to, less than or greater than 1, respectively. This kind of transition behavior for MSD is not 
unprecedented. In polymer physics, for example, it has been observed that the MSD of a single 
bead in a thermal fluctuating bead-spring chain will undergo MSD transitions |6CH 138) . 

In single particle tracking techniques, the local MSD is typically estimated by means of a 
heuristic regression-based method, which we will call pathwise MSD. Suppose that a microrheo- 
logical experiment generates a tracer bead sample path with observations X(Aj), j = 0, 1, ...,N. 
Let fj-2(t) ■= EX 2 {f). Under (fT7|) . assuming the increments are stationary, for h = 1,2, ...,n and 
n « N, one hopes that 



E [X 2 {t)] « at a , a > 0, tel. 



(17) 



1 



N-h 



/2 2 (A/i) : 



N-h + 1 



(X(A(j + h))-X(Aj)) 2 ^E[X 2 (Ah)] , if N is large. (18) 



j=0 



One then considers the linear regression specified by 

log(/2 2 (A/i)) = log(cr) + alog(A/i) +e h , h= l,...,n, 



(19) 
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Figure 4: At left, a plot of the pathwise MSD for 100 Brownian paths with a diffusivity that 
corresponds to a 200 nm radius spherical particle in water at 20° C. At right is the same pathwise 
MSD displayed on a log-log plot, which is prominent in the literature. Using a log-log plot hides 
the considerable variation in the empirical estimate for the viscosity of the fluid. Furthermore, a 
substantial number of paths will appear subdiffusive when subjected to algorithmic analysis. (See 
also Figure [5j) 



and applies ordinary least squares to obtain the estimator (log(cr), a). 

However, this procedure is subject to several shortcomings. First, the problem of delineation 
of regimes is inherently ill- formed, so the methods are naturally ad hoc. One must always "eyeball 
it" to decide when the MSD curve looks approximately straight on a log-log plot, and typically 
the transitions between regimes only become evident when an entire ensemble of pathwise MSDs 
are superimposed. 

As mentioned in the last section, determining whether particle paths are asymptotically dif- 
fusive or subdiffusive is a fundamental dichotomy in interpreting experimental results. However, 
data for large time lags is scarce and when available it may be overwhelmed by spurious drift in 
the fluid medium. Even after correcting for these physical constraints, the error between pathwise 
MSD and the true diffusivity of the particle is itself subject to high variability and induces de- 
pendencies as an artifact of the technique. The problem arises because the error terms in (1191) 
need to be at least approximately independent and identically distributed for standard regression 
properties to apply. However, this assumption is violated in the cases of interest. Indeed, as the 
lag becomes larger, the variance of /22(A/t) tends to increase. Also, neighboring (and even distant) 
estimators ^(A/t) are formed from the same observations, thus implying that they are dependent 
even in the diffusive case. Moreover, the convergence speed of the estimators /^(A/i), as a func- 
tion of N can be slow due to the dependence structure of the sequence (X(A(j + h) — X(Aj)) 2 , 
j = 0, N — h, in both diffusive and subdiffusive situations. 

The inaccuracy of the pathwise MSD can be at least partially compensated by the use of large 
amounts of data. However, the number of observations of laboratory generated tracer bead paths 
is usually limited to no more than 5,000. The left plot in Figure [7] illustrates the performance of 
the pathwise MSD based on Monte Carlo runs of Brownian motion sample paths (a = 1). Even 
though the increments of the underlying process are actually independent, the histogram for the 
distribution of the a as estimated by the pathwise MSD is concentrated in the very wide range 
[0.4,1.4]. 

Another serious drawback to the pathwise MSD based methods is that they lack an obvious 
way to quantify sampling error. For instance, a natural technique would be to perform parametric 
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bootstrap. One assumes that the estimated MSD is the true value and then generates correspond- 
ing paths. However, there are multiple memory kernels giving rise to the same theoretical MSD, 
and thus it is impossible to capture all the stochastic behavior of the underlying stochastic process 
based on an estimate of a alone. Nevertheless, the ability to quantify estimation error is especially 
important when trying to detect heterogeneities across multiple paths in a sample. The inabil- 
ity to quantify sampling variation has been a major impediment to the development of testing 
procedures for microrheological inferential techniques. 

3.2 The Diffusive Case 

Assuming one can infer that a process is asymptotically diffusive, the focus of statistical investi- 
gation should be the transient anomalous behavior. Because the rheological statistic of interest - 
the complex viscosity rf - is essentially the Fourier transform of the memory kernel T, it suffices to 
directly pursue the memory kernel in the time domain when possible. Fricks et al |25j developed 
such a method in the case where the memory kernel is assumed to be a Prony series; recall equa- 
tion (|lip . This is the class of memory kernels associated with well-studied Maxwell fluids [24j . 
and therefore any microrheological methodology should be able to reproduce characterizations 
from classical rheology. 

The basic idea of the method is to use a higher dimensional linear stochastic differential 
equation representation of the generalized Langevin equation, a fact that was previously developed 
in several forms \27\ [52] . Since the resulting process is a finite-dimensional linear ltd stochastic 
differential equation, the conditional mean and variance can be calculated for the sampled data 
permitting the use of standard Kalman methods to derive a likelihood function on even coarsely 
sampled data. 

One of the benefits of this Kalman method is that standard likelihood theory applies giving 
asymptotic estimates of standard errors of the parameters. Given the likelihood function, a 
Bayesian method could also be devised allowing for similar quantification of the variability of the 
estimation scheme. The ability to quantify estimation error is especially important when trying 
to detect heterogeneities across multiple paths in a sample. 

Kou and Xie |37| used a version of the GLE to study the subdiffusivity of a particle attached to 
a protein; however, their methods could be adapted to the diffusive case. In particular, they used 
a form of the GLE with a quadratic potential [50j . By using Kou's representation of the covariance 
of the resulting equation, they were able to use a method of moments type estimator to fit the 
Hurst parameter to experimental data [36J. One could attempt to modify this method using 
a similar representation of the covariance of the increment process for the differenced position 
data to apply Bayesian or maximum likelihood methods to the resulting Gaussian process. Note, 
however, that this strategy could be computationally very intensive and its success could largely 
depend on the form of the memory kernel selected. Both of these are examples of methods to 
fit microrheological data in a modern statistical context while maintaining fidelity with the well- 
established physical models, though both use simplifying assumptions to yield a more tractable 
model. 

3.3 Detecting Subdiffusion 

Exploratory investigation of experimental particle path data on any hydrogel, and mucus in par- 
ticular, indicates that the qualitative behavior of the submerged particles depends strongly on size 
and surface chemistry. Even among apparently identical particles, significant variation is seen, 
presumably because mucus has a highly heterogeneous mesh size distribution that interacts dif- 
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Series diff(paths[, 1]) 



Figure 5: The ACF of a single path of 1800 observations taken of a particle in human lung mucus 
with 2% mucin from the Hill lab, representative of a healthy patient. The sample was prepared 
according to the procedure presented in [38] . The ACF lacks a clear signal of dependency between 
increments. Figure [6] will demonstrate that this data exhibits subdiffusivity, which cannot be 
detected here. 



ferently with particles of different radii. Typically these properties are reported in terms of highly 
variable pathwise MSD curves, and in terms of labeling subpopulations of a particle ensemble as 
diffusive, hindered or immobile. This pathwise classification, which has important consequences 
in applications to drug delivery, poses a very difficult statistical question: Is it possible to success- 
fully label a given path as being from a diffusive or subdiffusive ensemble, and how much data is 
necessary to accomplish this task reliably? 

For the purpose of capturing subdiffusivity in the data, it is worthwhile shifting the perspective 
from the time domain, typical of the currently employed heuristic methods (see (|19p ). to the 
Fourier (frequency) domain. The latter is often superior for the characterization of the long term 
memory of stochastic systems. Given sufficiently fine data, i.e. high resolution in space and time, 
it should be theoretically possible to extract this memory information from the increment process 
of the position data. However, the fundamental limitations of the camera frame rate prevent 
the collection of adequately fine data for characterization using the autocovariance function. For 
example, in Figure [5] one sees virtually no signal in the sample autocorrelation of the increment 
process. 

The basis for spectral domain inference for long-term subdiffusivity is the connection between 
the spectral properties of the velocity process V and the population quantity of primary interest, 
namely G*(u). More specifically, one can relate the large time MSD exponent a in (|17j) to the 
behavior of the spectral density of the velocity process at the origin, and in turn relate this to the 
behavior of the complex shear modulus G*(uj) near the origin. 

We assume that the spectral density p(u)) obeys a positive power law near the origin, i.e., 

-f^ = c>o, o<d< 1 - 

uj^o \oj\ 2d 2 



lim f-^ = C > 0, < d < -. (20) 



As is done in the long range dependence literature, we use the fractional parametrization d out 
of notational convenience, and also because it naturally connects with the Hurst parameter when 
applicable (see expressions (|29p and (|30p below). Under a mild assumption, one can show that 

E\X 2 (t)] , 
c < ^_;/ J < C, for large t, (21) 
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for two constants c, C > (see the Appendix). In particular, if there exists < a < 1 such that 
(|17p holds for large t, then 

a = 1 - 2d, (22) 

which establishes the connection between the spectral density and subdiffusivity parameters. In 
particular, the velocity particle of a GLE driven by fBm is an example of a process whose spectral 
density satisfies (|20j) (see Section B~2j) . 

We still need to explicitly connect the subdiffusive exponent and by extension, the parameter 
d to the complex shear modulus G* . Recalling that p(uj) oc l/r]*(uj) = iui/G*{ui) and putting this 
together with (|22p . we have the simple relationship between the asymptotic MSD exponent a and 
the observed behavior of the complex shear modulus G*. Thus, when there is a clear asymptotic 
MSD exponent < a < 1 we have 

Um = C > o. (23) 

For perspective, we draw attention to Figure El which contains storage and loss moduli measure- 
ments from the original MW [47] work. The subdiffusive exponent exactly corresponds to the 
behavior of the storage and loss moduli as the frequency approaches zero. 

In reality, experimental data sets are made up of discrete measurements of the position 
process X(t), not of the velocity V{t). We should therefore consider the increment process 
Yj = VX(Aj) = X(Aj) — X(A(j — 1)) associated with the position process X(t). To simplify 
the following expressions, we will assume that A is one; this is equivalent to reparameterizing the 
problem into time units of the sampling increment. It turns out that, under assumption (|20|) and 
an additional mild condition, Yj inherits the fractional behavior of V(t). In other words, let py 
be the spectral density of the discrete increment process Y. Then 

lira = C > (24) 

for some constant C; see the Appendix. The property ()24|) implies that the discrete increment 
process Y is antipersistent and that it tends to display negative correlation across short time lags. 

The expressions (|22p and (|24p imply that a natural first approach to the problem of charac- 
terizing long term subdiffusivity is via spectral semiparametric methods. Indeed, the estimation 
of the parameter d does not necessarily involve the estimation of the entire function p(uj), or the 
specification of a full parametric model for V(t). Instead, one only needs to look at the behav- 
ior of p(u) close to the origin, a common idea from the literature on long memory processes. 
The fractional property (|24h opens the possibility of Fourier or wavelet-domain methods such as 
the Local Whittle estimator (Kiinsch [3D], Robinson [58]), spectral regression (Robinson (5DJ) or 
wavelet regression in the fashion of Veitch and Abry [68] , Moulines et al. [53] , Bardet [5j |6j (see 
also Beran [8] and Palma [55] for surveys). In particular, the Local Whittle is a well-studied 
semi-parametric estimator defined as 

/ 1 -s~\ I n (u)j) \ 1 

argmm d6e log — > — + 2d— > log(wj), m < n, 25 
\ i=i 3 ) j=i 

where = [di,^], < d\ < tfe < 1- In (|25p . I n {') is the periodogram, obtained from the discrete 
Fourier transform of a discrete sample path Y\ , Y-i , .. . , Y n as 
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Figure 6: Estimated distribution of the Local Whittle estimator over various sample paths (data 
from the David Hill laboratory, UNC). The red line is a Gaussian curve centered at the local 
Whittle estimate from the pooled data. The black line is a kernel density estimate summarizing 
the local Whittle estimates across multiple paths. Left plot: mucus. Right plot: water. 



The Local Whittle makes mild assumptions on the underlying process, i.e., that it is linear with 
innovations that are martingale differences satisfying a uniform integrability condition (Robinson 
[58j , pp. 1633-1634). There is the additional tuning parameter m, which indicates the number 
of frequencies used in the estimation. This tuning parameter is necessary because the spectral 
density of the process is semi-parametrically assumed to follow a power law (|20p close to zero. 
So, the use of frequencies which are too distant from the origin could then introduce undesired 
information from other parts of the spectrum of the process. The Local Whittle estimator is 
asymptotically normal, and combined with expression (|22p this implies that 

y/m(a-a) ~ N(0,1), (26) 

where m satisfies 1/m+m/n — > oo as the sample size n goes to infinity, which allows for asymptotic 
hypothesis testing. 

The right plot in Figure [7] displays the simulated distributions of the pathwise MSD and the 
Local Whittle estimator for Brownian motion. Both are centered at about the right value a = 1, 
but the latter exhibits a substantially smaller dispersion. Moreover, preliminary results with 
experimentally collected data indicate a substantial discrepancy in the distribution of the Local 
Whittle estimator between the cases of mucus and water (see Figure [6]). This stands in stark 
contrast with the results obtained from the sample autocorrelation function, displayed in Figure 

El 

3.4 Challenges in Statistical Inference. 

There are a number of outstanding challenges in statistical inference for microrheology. One 
important factor in the analysis and modeling of microrheological data is the relatively complex 
nature of the experiments. In general, the tracer beads cannot be observed precisely. The primary 
data resulting from a microrheological experiment is recorded in the form of a digital video. 
The latter is preprocessed using off-the-shelf software such as Spot Tracker and the NIH-funded 
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Figure 7: Estimation of the subdiffusivity parameter a = 1 for Brownian motion (path length 
5,000; only the first 1000 lags are used). Left plot: pathwise MSD. Right plot: Local Whittle. 
The red line in the Local Whittle plot represents the theoretical sampling distribution centered 
at the pooled estimate for a. 



ImageJ to identify and track the location of each bead in the complex fluid sample. Consequently, 
separating out the errors generated by image analysis is a considerable statistical task in its own 
right. 

Another important challenge is drift removal. In a microrheological experiment, a drift can 
be induced by multiple factors including a possible heat gradient from the microscope lighting. It 
has yet to be determined whether the drifts felt by particles are different, the same, or possibly 
random. When developing methods to estimate either the memory kernel or mean squared dis- 
placement, one should at least consider the effects of such drift. By using standard exploratory 
tools, the authors were able to identify numerous anomalies in the data emerging from David 
Hill's laboratory such as non-linear drifts and mis-tracked beads which allowed the laboratory to 
refine their experimental data collection. It has been an ongoing challenge to find data sets which 
are sufficiently unaffected by such experimental issues to properly apply some of the theoretical 
tools that are being developed. 

There is also the issue of the detection of heterogeneity, or non-stationarity, in the (sub)diffusive 
behavior of beads within experiments. Biofluids such as mucus and cytoplasm have a highly com- 
plex chemical composition with at least 250 different chemical component which may contribute 
to such heterogeneity. Do different beads behave differently within the same physical sample? Is 
there spatial variation in the diffusive behavior for the same bead? The former could indicate 
variations in the way beads interact with the material, for example due to surface chemistry; the 
latter may indicate heterogeneity of the material under study, which may be difficult to discern 
and quantify visually. An outstanding issue is the modeling and inference for large sets of tracer 
beads taken within the same experiment (i.e., physical sample or locus), or across different paths. 
One could consider adapting "longitudinal" methods from the time series literature to quantify 
heterogeneities across paths, such as those found in Shumway and Stoffer [63] and Durbin and 
Koopman [22] in combination with the Kalman approach of Pricks et al [25 j . 

One important question is whether the GLE is the right model or even the right type of model. 
Some data has suggested that particles bind to the polymer network and then break free to then 
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bind again after diffusing some distance [62]. These type of models go by several names depending 
on the field: regime switching in economics, Markov modulated processes in applied probability, 
and a flashing ratchet in some sub-disciplines of physics. There are a number of methods to 
handle this type of switching behavior statistically [63J, [TUl E3] • However, one possibility is that 
particles exhibit subdiffusion when free from the polymer network, and this would present new 
challenges. 

Another challenge is to test and perhaps drop the implicit assumption in the literature that 
the two spatial dimensions of the bead's location (X%(t),X2(t)) are independent. One natural 
way to do this is to develop a new GLE-like model driven by the multivariate process operator 
fBm (e.g., Didier and Pipiras [19] ) and test it based on recently-developed multivariate estimators 
(e.g., Abry and Didier [I]). 

4 Stochastic Simulation. 

4.1 Methods of Simulation. 

There are various assumptions under which the GLE yields stationary solutions, either for the 
position process X(t) or the velocity process V(t). Stationarity allows for the application of a 
wide range of techniques designed for such processes under the assumption of Gaussianity; this is 
especially convenient since the system is in general non-Markovian. 

As described in [7], there are various ways to simulate a Gaussian process: (SI) if the co- 
variance is known, by means of a sequence of variables whose covariance is the given one; (52) 
based on a representation of the process (such as a time or spectral domain stochastic integral) , 
possibly after some simplification (e.g., truncation) of this representation; (5*3) through a discrete 
sequence whose weak limit is the target process. 

With regard to (S3), if a natural discretization scheme is available, simulating a continuous- 
time process is potentially straightforward. A GLE is a kind of SDE. In the case of ltd diffusions 
driven by Brownian motion, Euler-type schemes, which rely on discrete approximations at equidis- 
tant grid points, can be used for simulation. For instance, see [3] or [35]. In general, for SDEs 
which are driven by other types of Gaussian noise, the theory is not as well-developed; on the 
pathwise approximation to SDEs driven by fBm, for instance, see [54] and references therein. 
However, it is well-known that an equally-spaced sampling of the OU process produces an AR(1) 
sequence. Thus, a fast and exact simulation procedure is to generate an AR(1) vector with ap- 
propriate parameter values. In [25], the multivariate equivalent of this fact is used to simulate the 
solution to a generalized Langevin equation when the memory kernel is a Prony series. A related 
example can be found in [49], where the authors show that the same Prony series GLE can be 
represented as a Brownian motion plus the sum of OU processes when taking a zero mass limit. 
In general, though, the correlation structure of discretized continuous time processes, such as the 
GLE, is complex and not amenable to a simulation scheme based on a simple loop, as with the 
OU process. 

The availability for the GLE of integral representations or covariance functions of the form (1131) 
and (|14p . respectively, indicates that it is natural to pursue simulation methods that fall under 
(SI) or (S2). We now describe two such simulation methods which assume known covariance 
functions and that will be of interest in the simulation of GLEs. Later on, in Subsection 14.31 we 
will also describe a wavelet-based method that falls under (S3), but which builds upon integral 
representations. 

So, let X be a target zero mean Gaussian stationary process with a given, known covariance 
matrix S over a finite set of time points N. One classical and simple way to simulate X consists 
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of performing a Cholesky decomposition of £ = LL*, where L is a lower triangular matrix, and 
generate a vector Z of i.i.d. standard Gaussian variables. Then X = LZ, as desired. Cholesky- 
based simulation is exact up to computational error and can be implemented recursively (see [3], 
pp. 311-312). However, it is slow in terms of computational complexity: 0(A 3 ), where N denote 
the length of the resulting stochastic vector. 

Another popular method is the Circulant Matrix Embedding (CME; see [15] , [69], [20]; see 
|33j, [Sj and [3] for reviews of this and other methods). The algorithm draws upon embedding 
the covariance matrix in a non-negative definite circulant matrix of size M > 2{N — 1). This is 
computationally convenient, since the diagonalization of circulant matrices can be carried out by 
means of the Fast Fourier Transform (FFT), which has complexity 0(./Vlog(-/V)). Like Cholesky- 
based simulation, CME is exact. For the reader's convenience, we summarize the CME. Let J- be 
the FFT. The simulation procedure is as follows: 

(Ci): choose M = 2 P > 2(N - 1); 

(C 2 ): form the vector A := (A(0), A(l), A(M/2 - 1), A(M/2), X(M/2 - 1), A(l)), where \{h) 
is covariance function evaluated at the time lag h; 

(C3): compute A := ^(A). Check whether the entries of A are all non- negative. If so, proceed to 
the next step. If not, choose a larger p; 

(C 4 ): generate Z := X + iY where X, Y 2V(0, /(l/v 7 ^)); 
(C 5 ): define U, U k := sJT k Z k) k = 1, ...,M; 
(C 6 ): compute X := ^J 7 " 1 ([/))• 

The vector X obtained is Gaussian and has the desired covariance structure. Of course, the non- 
negativity condition in (C3) is essential. Although no general theoretical results are available, it 
has been justified for many models of interest. For instance, it is shown in [2D] that a sufficient 
condition is that the sequence of covariances A(0), A(l), X(N — 1) be non-negative, decreasing 
and convex (see also [3], [56], [13]). One limitation of the CME is that it is not iterative in the 
sense that it requires the simulation horizon to be predefined. 



4.2 The fractional GLE: correlation structures via integral representations 

In order to provide a study of simulation of GLE-based anomalous diffusion, we focus on the fBm- 
driven GLE (fGLE), since its correlation structure has been well-studied. In this section, based 
on [36], we take a closer look at the correlation structure of the fGLE via integral representations 
(see also Section [2T3|) . 

In that paper, the author analyzed the fGLE in three situations: when it models a free particle, 
when the particle is subject to the harmonic potential U (x) = ^mifjx 2 , and when the acceleration 
term mdV{t)/dt is negligible, corresponding to the overdamped condition in physics, where the 
friction is quite large. Because it is the most relevant for passive microrheology, we will focus on 
the free particle system, i.e., 



mdV(t) = -j( / V{u)T H {t - u)du)dt + y/ '27 k B TdB H (t). (27) 
J — 00 
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Figure 8: fGLE, correlation structure up to a multiplicative constant (d = 0.25, j = 2, m = 1) 
Left plot: spectral density. Right plot: covariance function. 



Kou [36] gives the solution V(-) to (|27p as a pathwise defined Riemann-Stieltjes integral. In the 
spectral domain, one can develop such solution into the integral representation 



v(t) c V^T , ^ 



1 



1 



7k(w)|w| 



-2d 



imuj 



\to\~ d dB(uj), 0<d<-, 



(28) 



where 

d:=H- 1/2, (29) 

7 is the friction constant, ks is the Boltzmann constant, T is the underlying temperature, c(d) 
is a constant, B(u) := B\{u) + ii?2( w ) for the real-valued Brownian motions B\, B<i satisfying 

dB(-uj) = dB{uj) a.s. and 

k(uj) = T(2H + l)(sin(F7r) - isign(w) cos( J fTvr)). 
We can rewrite the spectral filter of V with respect to the complex Brownian measure dB(-) as 
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/?:=! + 2d, 



(30) 



for appropriate constants i/q, v\, v<i- Expression (|30p is the basis for the simulation of fGLE in 
Subsections 14.31 and 14.41 In particular, expression ([28]) implies that the fGLE is antipersistent 
(see (|20p ). Figure [8] depicts its correlation structure in the spectral and time domains. 

4.3 Wavelet-based simulation 

In this section, we describe a new, approximate wavelet-based simulation method for the fGLE. 
The discussion is based on [17] and [18], where the mathematical details are developed. 

The simulation method goes under (53), since it generates sequences of converging discrete 
time processes. This method builds upon approximate discretizations generated according to 
analytical and computational convenience. Moreover, it has the following features: it is com- 
putationally fast, potentially reaching complexity O(N), since it is based on a Fast Wavelet 
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Transform-like algorithm; it provides approximations that converge uniformly over compact in- 
tervals almost surely; the convergence speed is exponentially fast and depends on the sample path 
smoothness of the limiting process. Moreover, it is iterative both intensively and extensively. In 
other words, a generated approximation at scale J over a compact interval [0, T], T £ N, can 
be used to generate a finer approximation at scale J + 1 over [0, T] or some expanded interval 
[0,T + r], r £ N. 

The simulation method involves an appropriate sequence of filters gj, j = 0,1,..., J, where 
J is the finest approximation scale chosen, and an associated sequence of low- and high-pass 
wavelet filters Uj and Vj, respectively. The method then amounts to recursively generating an 
induced sequence of discrete approximations Vj, at scale j 6 N, of the continuous time velocity 

process V, where Vj : k = Yl^=-oo 9j,nCk-n, and '~ ' WN(0, 1). The choice of the sequence of 
discretization filters gj is the key requirement. Heuristically, it should be such that 

Gj(2~ j u) := ft^f* ~ weR, (31) 

where gj{oo) is periodically extended to R. Intuitively, the discretization filter approximates, up 
to a scaling factor, the continuous time filter as the scale becomes finer and finer. 

So, let f2 be the upsampling by factor 2 operator, i.e., it inserts a zero between any two entries 
of a vector. The wavelet filters Uj and Vj are defined in the Fourier domain as 

%M := ^-TTTT Vj(uj) :=g j+1 {uj)v(uj), (32) 
gj{zuj) 

where u, v are the conjugate mirror filters (CMF) of the underlying wavelet multiresolution 
analysis (see chapter 7). The simulation procedure can be briefly described as the following 
Fast Wavelet Transform-like algorithm. 

(Wi): at scale/step 0, generate via an exact method (e.g., CME) one first approximation sequence 
Vb; 

(W2) : at scale/step j € N, and given an approximation sequence Vj, obtain the next discrete 
approximation V^+i at scale/step j + 1 via the relation 

V j+1 = Uj* | 2 Vj + Vj* |2 £j- 

The algorithm works because at step j + 1, the FWT annihilates the correlation structure 
(kernel) gj of the approximation { V^} at scale j, and replaces it with a new, pre-chosen correlation 
structure ffj+i- This is a consequence of (p?2"j) and of the properties of the CMFs u and v (see also 
[57]). 

As a simulation procedure, the resulting wavelet coefficients converge appropriately to the 
stochastic signal, i.e., 2 J / 2 Vj[ 2 -'t] ~~ ^ V(t), where [x] denotes the integer part of x E R. Such 
property is a consequence of (|3ip and of the properties of the scaling function <j). The choice of 
the function Gj(-) plays a key role. If V is an OU process, then equally sampled points make up 
an AR(1) sequence. Therefore, for simulation purposes, 

(I _ e - A2_J e -^)- 1 

G ^ = k i ■ U ( 33 ) 

(A + 1 

is a natural choice, since (1 — e _A2 ' 'e^)^ 1 is a spectral filter of the associated AR(1) sequence. 
Moreover, the heuristic relation ([3T|) is indeed satisfied. 
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However, exact and simple discretization schemes such as (|33|) are not generally available. 
As a rule, the discretization of a continuous time process leads to a substantially more intricate 
expression for the spectral density. This is true, for instance, for fBm and the fGLE. In particular, 
for the latter, even though expressions for the covariance structure of the fGLE are available, they 
do not appear in closed form, which calls for numerical methods. 

So, a natural question is whether one could construct converging discretizations by means of 
a simple method that applies to a wider class of stochastic processes, in particular, the fGLE. 
Indeed, this can be done by developing non-causal discretization filters gj in three elementary 
steps: 

(t.l) extend the truncated function ^(w)lr w n \ periodically to R (i.e., gj stems directly from g); 

(t.2) modify the resulting function with rescaling terms (e.g., 2~- J ) so that relation (|3ip holds; 

(t.3) smooth the resulting function at — tt, tt so to speed up the time domain decay of the filter 
in theory and computational practice. 

The method described in steps (t.l) — (t.3) produces a sequence of discrete time processes Vj = 
gj * e which is approximate in two senses. First, because of numerical error including truncation, 
which is inevitable in computational practice with convolution-based procedures. Second, because 
gj is picked for analytical and computational convenience, and thus Vj is not in general an exact 
discretization of V. However, one can show that these approximate discretizations still converge 
exponentially fast to the correct limiting process. 

By applying (t.l) — (t.3) to (f30j) . the proposed filter for the free-particle fGLE is 

9j(A =9j,iyA UJ )9jA UJ )i ( 34 ) 

where 

?,, d H = 2^exp(^(^) 2 ) u* d :=irVd, ue[-tt,tt). (35) 



TT 



(A 



i 



g hU ^) := exp ^ j - - - ~^^ 2 , - € [-*,*). (36) 

In (|35p and (|36p . the exponential terms smooth the original truncated filters. Figure[9l right plot, 
shows the designed high-pass filter in the time domain. The oscillations are almost invisible for 
relatively low lag values, and the actual decay is fast. Moreover, <7j(2~ J cj) ~ g(u) for large j, as 
expected. For the sake of comparison, Figure Figure El left plot, displays the time domain filter 
obtained without the smoothing terms. 



4.4 Evaluating simulation methods 

Most simulation techniques are supported by theorems that establish some sort of convergence, 
equality in law, and so on. However, the finite sample performance, or accuracy, can in principle 
be disparate across methods. One idea is to use estimators in order to compare the simulation 
methods. Nevertheless, only the asymptotics of estimators are available in most cases, and the 
finite sample performance of estimators, e.g., bias, is studied based on simulation, which creates 
a circularity. 

In view of this, we can look at measures of performance of the methods relative to one another. 
Since Cholesky-based simulation is a simple and exact procedure, we choose it to provide the 
baseline for the CME and the wavelet-based method (see also [7], [E]). A two-sample t statistic 
is used to assess the difference between the values of the estimator when generated by two of the 
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Figure 9: Left: Non-smoothed high-pass filter, time domain. Right: fGLE, high-pass filter at 
j = 4 (7 = 2, m = 1, d = 0.25), 4 zero moments, time domain. 
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methods. We evaluate the quality of the simulation based on the Local Whittle estimation of the 
parameter d. This is of special interest in our framework of subdiffusion, since the Local Whittle 
is a good estimator for the subdiffusivity parameter, as discussed in Section [3l 

Since the direct discrete sampling of the fGLE {V(i)}teN does not generate an antipersistent 
process, for the purpose of comparison we opted for simulating the Cholesky and CME sequences 
based on the covariance function of the increment process Y% = VX(i) = X(t) — X{t — 1), i.e., 

E[Y t Y t+h ]=c*(d) [ e 1 ' 

JR 

where T = 2 9 , and c*(d) > is a constant. In turn, the wavelet simulation consisted of generating 

the fGLE, approximating X(t) ~ X^L=/ T ^' Vj,k{T/2 J ) and then sampling over 1,...,T. The 
wavelet filters are truncated at lag \L\ = 80, which ensures that at the point of truncation the 
filters Uj and vj, j = 1, 13, exhibit terms of magnitude below 10 -5 . The autocovariances 
(|37p and the integrals of the initialization sequence (at j = 0) of the wavelet method were both 
calculated using quadrature-based numerical integration. 

The results are shown in TablesQ]for d = 0.10, 0.25, 0.45. In all cases, the relative performance 
of the wavelet method seems to be fairly insensitive to the final approximation scale J. This 
indicates that it is not necessary to use a very fine scale in order to attain higher accuracy 
levels. On the other hand, since the wavelet method is the only approximate one among the three 
considered, it is a priori expected to be more prone to affect the bias of the estimator. Indeed, 
the results seem to indicate that Cholesky and CME-based simulation tend to be slightly closer 
to each other, as measured by the estimator bias. Nevertheless, all three methods seem to be 
close enough, i.e., within two standard deviations. 



i 



Un -I- Vl \uj\P A- ZA> 



\co\ 2d dio, 



t = l,...,T, (37) 
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Table 1: fGLE: Local Whittle estimation of the parameter d (7 = 2, mass = 1), comparison across 
values of J, AX(n), time series length 2 9 

Uj, Vj cut off at lag \L\ = 80, init. CME time series length 2 10 
d = 0.10 



method 






d 


s 


JV 


\t\ statistic 


wavelet (CME at j = 0, J 


= 6) 




0.12044083 


0.09245007 


5000 


1.05431590 


wavelet (CME at j = 0, J 


= 8) 




0.12241433 


0.09102200 


5000 


0.02227624 


wavelet (CME at j = 0, J : 


= 10) 




0.11932462 


0.09357999 


5000 


1.65282860 


CME 






0.12221790 


0.09220211 


5000 


0.08515547 


Cholesky 






0.12237381 


0.09088329 


5000 




d = 0.25 


method 






d 


s 


JV 


\t\ statistic 


wavelet (CME at j = 0, J 


= 6) 




0.27857299 


0.09217255 


5000 


0.25582876 


wavelet (CME at j = 0, J 


= 8) 




0.27998323 


0.09268102 


5000 


0.50777309 


wavelet (CME at j = 0, J : 


= 10) 




0.28029380 


0.09627634 


5000 


0.66273074 


CME 






0.27818405 


0.09112314 


5000 


0.46947509 


Cholesky 






0.27904459 


0.09144735 


5000 




d = 0.45 


method 






d 


s 


JV 


\t\ statistic 


AWD (patched filter, CME at j - 


= 0, J 


= 6) 


0.43156246 


0.10538368 


5000 


0.78943448 


AWD (patched filter, CME at j = 


= 0, J 


= 8) 


0.43199573 


0.10660877 


5000 


0.58088823 


AWD (patched filter, CME at j = 


-- 0, J - 


= 10) 


0.43166818 


0.10678505 


5000 


0.73448769 


CME 






0.43338219 


0.10403330 


5000 


0.07274950 


Cholesky 






0.43322954 


0.10579033 


5000 





4.5 Challenges in Simulation. 

Simulation procedures are important to assist inferential efforts in multiple contexts in microrhe- 
ological research. Ideally, one is interested in fitting models based on memory kernels taken from 
a broad class of functions. It is still an open question whether CME and wavelet methods are 
well-suited in such a general context. This is of interest, for instance, in Monte Carlo-based 
quantification of sampling error. 

Once one has fit a model for the diffusive properties of a material, there are a range of 
applications related to exploring the physical ramifications of the model. In particular, there has 
been strong interest in characterizing the properties of human lung mucus to better understand 
the harmful effects of cystic fibrosis on lung health. One application of Monte Carlo methods 
is the study of first passage times to investigate drug delivery vehicles. When a drug delivered 
via small particles is introduced into the lungs through such devices as inhalers, the particles 
will diffuse on the surface of the lung through mucus. By studying the first passage times of 
the particles, one gains a better understanding of the distribution of the particles throughout the 
internal lung surface. A related challenge is to simulate diffusing particles in a properly modeled 
mucus layer, but with complex boundaries determined by the internal lung surface. The particles 
will be reflected at this possibly non-smooth interface, a complex effect that needs to be taken 
into account in simulation efforts. 
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Two point microrheology expands the challenge of simulation to multiple coupled beads. In 
|32j, the focus is on the hydrodyamic interaction of two beads in close proximity. However, this 
effect does not have to be limited to only two. One of the reasons for using two is to simplify 
analytical results, and simulation would allow one to move beyond this limitation. A related issue 
is to devise simulation methods for the full generalized Langevin equations including a possibly 
non-quadratic potential ([9]). The existence of solutions of the stochastic differential equation 
becomes much less clear; moreover, the typical methods for devising numerical schemes typically 
relies on an Ito-type formula, which may not be readily available. 

5 Appendix 

The relation between the behavior of the spectral density at the origin and the decay of the auto- 
correlation function of a stationary process outside certain classes of parametric models is subject 
to specific assumptions (see for the antipersistent case, see [U]). However, the related issue of 
the asymptotic behavior of the MSD can be addressed under general assumptions, as shown in the 
next proposition, used in Section 13.31 The statement is restricted to anti-persistence (condition 
([20]) ). but the same argument applies in the case of a singularity at the origin. 

Proposition: Let {V(t)}t>Q be a zero mean Gaussian stationary (velocity) process with spectral 
density p(uj) satisfying 

/•oo 

/ u? v log(l + u)p{uj)duj < oo (38) 

J 

for some v S (0, 1). Moreover, assume that p(co) satisfies $20\). Then $21\) holds. 
Proof: Let 

X(t) := [ V(s)ds (39) 
Jo 

where the integral (|39p is taken in the Lebesgue sense. From Cramer and Leadbetter [T3], pp. 
181-182, (I39p is well-defined in the sense that there exists a process rj(t) which is equivalent to 
V(t), and which satisfy a Holder condition of order v a.s. For notational simplicity, we will keep 
writing V(t). 

Let g(x) be a spectral filter for the velocity process V(t). Note that s, s' > 0, 

K[\V(s)V(s')\] < y/E[V\s)WE[V2(s>)], 

which is finite and constant, by stationarity. Thus, based on the spectral integral representa- 
tion of V(t), by applying Fubini's Theorem and to E[X(t)X(t% t,t' > 0, we obtain the 
representation 

c I ^W w)d g (w) (40) 

where B{oj) = B\{oj) + zi?2 ( w ) is a complex-valued Brownian motion such that dB{—uj) = dB(uj) 
almost surely. 

Condition (|20p implies that 

ci\u:\ 2d < p{u}) < c 2 \co\ 2d , \cj\<6, (41) 



26 



for some c%,C2,S > 0. Therefore, 



E [X 2 (t)] 



Atui 



ILO 



p(u))du) 



+ 

\w\<5 J\lo\>S 



Atui 



After a change of variables, 

2 



L 


e ituJ - 1 







uj\ 2d duo ~ t 



l-2d 



poo 


e ts - 1 


1 — oo 


is 



\ s \™ds = t 1 - 2 * E 



p{uj)duj. 



^V2-d(l) 



(42) 



where -Bi/2-d(") is a fBm with (Hurst) parameter 1/2 — d (see Section [2.3|) . Therefore, for large 
enough t, ([4*2]) can be bounded from above by 



c 2 t 1 ~ 2d E 



+ 



n>5 



2 



p(uj)duj, 



where the second term is finite because p(oj) G L 1 (M). On the other hand, for large enough t, 

^1/2— ^ or some constant c' x > 0. Therefore, 



2]) can be bounded from below by c^t 1 2d E 
2"Tj) holds with the constants c[, ci- □ 

The next proposition is used in Section 13.31 



Proposition: Let {V(t)}t>o be a stochastic process satisfying the assumptions of the previous 
Proposition. Additionally, assume that 

there exists a < 5 < 2tt such that p is bounded and continuous over (— 5, S) c . (43) 



Then holds. 

Proof: The discrete increments of (|4(jp have the integral representation 

g(oj)B(du), j G 



Yj = VX(j) 
For notational convenience, let 



lid 



1 - e 



h v (uj) = ^— ^ p(u). 

Then, by a standard calculation, the covariance of Yj is 

e iju ^2hy(u + 2Trk)du, j G Z. 
_7r fcez 

Therefore, the spectral density of Y can be expressed as 



ILC 



'p(oj) + |1 - e~ iul \ 2 Ku + lirk) 

fcGZ\{0} 



w + 2vr£:| 2 ' 



LO G — 7T, 7T 



We are interested in the behavior of the ratio Tjpr as w + . Note that, for small enough w > 0, 
w + 2-7r/c G (— 5, 5) c , k G Z\{0}. Therefore, by (|43[) and the Dominated Convergence Theorem, 

lim V p(co + 2iTk)- — — — ^ < 

fcez\{o} 

The claim then follows from (1201). □ 



oo. 
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